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Abstract 

An inverse problem of identifying inhomogeneity or crack in the workpiece made of nonlinear magnetic material is 
investigated. To recover the shape from the local measurements, a piecewise constant level set algorithm is proposed. 
By means of the Lagrangian multiplier method, we derive the first variation w.r.t the level set function and obtain 
the descent direction by the adjoint variable method. Numerical results show the robustness and effectiveness of our 
algorithm applied to reconstruct some complex shapes. 

Keywords:Nonlinear electromagnetism, adjoint variable method, shape reconstruction, piecewise constant level set 
algorithm. 
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In many applications, one needs to find out the flaws in materials nondestructively. Inspired by this, the non-destructive 
evaluation technique has attracted the eyes of many researchers during the last decade. As one kind of non-destructive 
evaluation technique, eddy current testing technique [T] has been used for flaw detection. In this paper, we intend to 
design an algorithm to identify the crack or inhomogeneities in the nonlinear magnetic material like steel from the local 
measurements of the magnetic induction. 

This inverse problem is a problem of shape reconstruction. Similar to general shape recovery problem, there is no 
information about the interface of the optimal shape as a prior, so we need to have a good mechanism to express the shape 
and track the evolution of the shape. The level set method was first originally proposed by Osher and Sethian in [5]. For 
this method, the interface between two adjacent domains is represented by the zero level set of a Lipschitz continuous 
function. Through the change of this function, this method can easily handle many types of shape and topological changes, 
such as merging, splitting and developing sharp corners. Due to these merits, it has been used in various areas, such as 
epitaxial growth [3], inverse problem, optimal design [3], image segmentation [S], structure topology optimization [5] and 
EIT problem [7|. For the sake of the numerical stability, the level set function is usually chosen as the signed distance 
function, but at most cases, the level set function after each iteration is not be a signed distance function and is usually 
re-initialized by solving an ordinary differential equation [5]. In [5], Cimrak et al. used the level set method to represent 
the shape of the inhomogeneity and evolve the shape by minimizing a functional during the iterative process. In jlOi lllj . 
Cimrak also used the level set method for the representation of the interface to solve some inverse problems in thermal 
imaging and the nonlinear ferromagnetic material. As to the initial value of the level set function, it was reported in 
[m [131 [2] that the level set method based only on the shape sensitivity may get stuck at shapes with fewer holes than 
the optimal geometry in some applications such as structure designs. When one wants to use the level set method to solve 
the practical problem, he can reduce the effects of the initial value of the level set function on the final results in the next 
two ways. The first way is to choose the shape with enough holes as the initial value. The second one is to introduce the 
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topological derivative into the level set method to let the shape create holes in the iterations [IH [151 HSl El IE] ■ To find a 
good initial value of the level set function, the researchers in [5] proposed the gradient-for-initial approach which is based 
on the idea that the domain which can drop the value of the cost functional should be the air gap. In that approach, the 
parameter in smeared-out Heaviside function should be set large enough. The large value of this parameter, however, may 
cause the oscillation phenomenon. So the parameter in the acquisition of the initial choice of the level set function and the 
evolving process of level set function should be set separately. 

Recently, piecewise constant level set method, which is a variant of level set method, was proposed by Lie, Lysaker 
and Tai in [5D1 HH [52]. To distinguish these two methods, we call the former one traditional level set method. Unlike 
the traditional level set method, the interface between two adjacent sub-domains is represented by the discontinuity of a 
piecewise constant level set function. Compared with traditional level set method, piecewise constant level set method has 
at least two advantages. One merit is that it can create many small holes automatically without the topological derivatives 
during the iterative process. Furthermore, it is verified by many numerical examples that the final result is independent 
of the initial value of the level set function in many numerical tests. And the other one is that the piecewise constant 
level set method need not to re-initialize the level set function periodically during the evolution process, thus, reduces the 
computational cost a lot. Since it was proposed, it has been applied in various fields such as image segmentation, elliptic 
inverse coeSicient identification, optimal shape design, electrical impedance tomography and positron emission tomography 
[19]. Lie, Lysaker and Tai took this method to solve the image segmentation in ^20, i2I, |22, ,23j and the elliptic inverse 
problem and interface motion problem in [211 [5S]. Wei and Wang used piecewise constant level set method to solve 
structural topology optimization in |26| . Zhu, Liu and Wu applied the piecewise constant level set method to solve a class 
of two-phase shape optimization problems in [27] . Zhang and Cheng proposed a boundary piecewise constant level set 
method to deal with the boundary control of eigenvalue optimization problems in [28]. In this paper, we intend to use the 
piecewise constant level set method to represent the shape and recover the exact shape of the crack or the inhomogeneities 
of nonlinear magnetic material. 

Based on the piecewise constant level set method, we propose a piecewise constant level set algorithm to recover the 
shape of the nonlinear magnetic material problem. We introduce the piecewise constant level set method and convert the 
constrained optimization problem into an unconstrained one by the Lagrangian multiplier method. In the numerical tests, 
our algorithm relies little on the initial guess of the level set function. Moreover, our algorithm can reconstruct shape 
accurately even when the noise level is high. 

The rest of this paper is organized as follows. In Section 2, we introduce the piecewise constant level method in brief. 
In Section 3, we introduce the direct problem and the inverse problem. In Section 4, we describe the deduction of the first 
variation w.r.t level set function for the objective functional in the unconstrained optimization problem in detail. In Section 
5, we present numerical results to show the effectiveness and robustness of piecewise constant level set algorithm. 



2 Piecewise constant level set method 

We first introduce the piecewise level set method in brief. Suppose that the domain is the union of some sub-domains 
Qi, i = 1, 2, • • • , m, i.e., 

m m 

^ = U ^» u U 

1=1 1=1 

where dfli is the boundary of sub-domain fli. 
If there exists a function defined as 

(j){x) = i, x^rii, where i 1, 2, • • • , m, (1) 

then the interface between two adjoint sub-domains can be identified by the discontinuity of the function and the 
characteristic function of the i-th sub-domain fli can be written as 

n (0 - j) /ON 

For the traditional level set method proposed in [2] , if is a level set function, then the expression of Xi contains the 
Heaviside function H(x), which is defined as below: 
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Since H{x) is not differentiable at 0, it is often replaced by smooth functions which contain some parameters in [25]. Some- 
times the parameters in the replaced smooth functions can cause difficulties in the result analysis, such as the convergence 
of the level set function </> in [3]. For the piecewise constant level set method, we do not encounter this kind of problems. 
If (j> is defined as in then it holds that 

K{(j>) = 0, in (4) 

where the function K is defined as 

K{<j>) = (0 - 1)(0 - 2) • • . - m) = - 0- (5) 

1=1 

It should be pointed out that if Eq. (jl]) holds, then every point a; G is in one and only one sub-domain. In other 
words, there is no vacuum and overlap between two different sub-domains. 

For any piecewise smooth function f{x) that coincides with fi{x) in fij, it can be written as 

m 
i=l 

where (j){x) is defined in ([T]). 



3 Direct problem and inverse problem 

Wc first introduce the quasi-linear partial differential equation which will be used in this paper. 
This equation is defined as 

r V ■ |VA|2)VA) = J, in O, 
\ yl = 0, on 50, 

where O C is a bounded domain with boundary, the symbol J denotes a suitable function defined in fl and the 
function tj : O x M — > K is defined as 

-(^'■^) = I ""ny^n (8) 

^ ' ' \ V2{s), X e n\D, 

where D C vi,V2 are two functions determined by the specific practical applications. 

Eq. ([7]) can model many industrial and physical applications. When both vi and V2 are constant functions, it is the 
elliptic inverse problem in |24| . When vi and V2 are both linear ones, it characterizes the electric impendence tomograph 
problem |29j . Since the material in this paper is nonlinear magnetic, we use the nonlinear model, which is more robust and 
able to reconstruct the shape of the crack or inhomogeneities even when the linear model is not [9]. In this model, V2 is a 
nonlinear function. 

Since we need to solve the direct problem in each iteration, we next present the direct problem in detail. The direct 
problem is to obtain the quantity A when the variable v, the domain D and the function J are given. In this paper, our 
goal is to design a numerical scheme to identify the crack or the inhomogeneity in the workpiece which is made of nonlinear 
magnetic material. In this case, the function J denotes the induced current density, the quantity A denotes the only nonzero 
component of the vector potential A which is perpendicular to the xy— plane. According to the physical knowledge, the 
magnetic induction B is 

B 

Since A is defined on a subregion of , we obtain 

After introducing two symbols fii and fi2 to denote the magnetic permeability of the air and the nonlinear magnetic 
material, we define the functions vi and V2 as the reciprocal of the two functions fii and /i2, respectively, i.e. vi ^ l/fii and 
V2 = ^2- As to /i2, due to the nonlinear magnetic property of the material, it depends on |VA|. Usually, the function 
/Z2(s) is monotonically increasing on the interval [0, Smax] and is monotonically decreasing on the interval [smax, +oo], where 
Smax is a number determined by the property of the magnetic nonlinear material. 



(-,--,of. 



|B| = |VA|, 



(9) 
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With the piecewise constant level set method, the function v{x) can be formulated as 

v{x)=vi{2-(b)+V2{<j}-l). (10) 
By ([TU]), the weak form of Eq. ([7]) can be written as 

/ [vi{\VA\^){2-(l)) + v2{\\7A\^){(f>~l)pA-V(pdx^ ( Jipdx, G VKo'^(17). (11) 
Jn Jn 

The existence and uniqueness of the solution to the Eq. ([TT]) . have been proved in [3D1 [HI]- To ensure the existence and 
uniqueness of the solution to Eq. (|lip . we adopted the assumptions in [S] which are listed as below: 

Al The function Vi is non-decreasing; 
A2 lims^ot'j(s) = v„iin > 0; 

A3 lims^oo Vi{s) =^ Vjnax > and define v{s) = Vmax for s = oo; 
A4 Vi is differentiablc with well-defined derivatives satisfying 

where z = 1, 2. 

In order to let the function /i2 satisfy the non-decreasing property in the assumption Al, we set the function J large 
enough. As to the solution to Eq. ([T]), we consider it as the solution to the following nonlinear operator equation 

G{A) = J, 

where G is an operator defined in the space WQ''^{Vt) and takes values in Wg^'^(ri). To solve this nonlinear operator equation, 
we first choose an initial guess Aq of A, and use the Newton-Raphson algorithm to update A, i.e., 

=A,- [DG{A,)]-\G{A,) - J), 

where i = 1, 2, • ■ • . 

Now we introduce the inverse problem. When vi and V2 are given, the inverse problem is to reconstruct the actual shape 
of D from the measurement data M of magnetic induction in a given subset F C il. Under the framework of the piecewise 
constant level set method, the problem is to find a piecewise constant function </> to approximate the exact level set function 
which can represent the actual composition of the workpiece (during the iterative process, the level set function (f> may be 
not a piecewise constant one). In order to distinguish these two functions, we denote by <f>exact the exact level set function 
here and afterward. For our work, F = fi. 

When the level set function and the data M on F are given, we use the following functional 

F,{4>)^\j^\'^A{c^)~M\^Ax, (12) 

to measure the misfit between (j> and (j>exact, where M = {—B2, Bi)'^ with Bi and B2 being the x-axis and y-axis component 
of B. Obviously, the smaller the value of Fi is, the more close the level set </> is to (pexact- Thus the inverse problem is 
converted into the following optimization problem 

minimi ((/«) subject to K{(t>) = 0, (13) 

where K{<t)) is defined in ([5]). 

4 Algorithm 

In this section, we will introduce our piecewise constant level set algorithm for this inverse problem. 
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Wc first give the deduction of tlie gradient of Fi. Introducing two symbols DFi and S^g to denote the gradient of Fi 
and the Gateaux derivative of g with respect to in the direction /i, we give the deduction of ShFi. According to the 
definition, it can be expressed as 

S,F, = lim + = / ^s,A . (VA - M}dx, (14) 

where VShA will be given below . 

Differentiating both sides of Eq. pip with respect to 0, we obtain the following sensitivity equation 

/ Sk[vix,\\/A\^)]\/A-Vipdx+ I v{x,\VAf)\/6hA-\/Lpdx=0. (15) 
Jn Jo. 

Since A is completely determined by the level set function (j) when J, vi and V2 are given, from (jlOp . the integral part 
of the first term in (fTSl) can be sorted as 



5h[v{x, |VAp)] =(t;2(|VAn - v,{\VA\^))h 

+ (2 - cj))5h[vi{\VA\^)] + {cjy- l)5h[v2{\VA\^)] 
={v2{\VA\'')-v^{\VA\^))h 

+ 2((2 - cj,)v[{\VA\'') + {cjy- l)v'2{\VA\''))VA ■ VSnA. 



(16) 



Substitute p6p into ()15p and simplify the equation, the sensitivity equation (jisp can be written as 

/ {v2{\yA\^) -v^{\yA\^))hyA-yLpdx 

Jn 

= - / 2((2 - <f))v[{\\/A\^) + {(f)- l)u2(|VAp))VA • VShAVA ■ Vipdx (17) 

- I v{x,\\7A\'^)\7ShA-Vipdx. 
Jn 

As to the computation of VfJ/jA, it could be deduced by the adjoint variable method which was used in many problems 
[32l l33l [34l [35l [36l [37]. Specifically, we take the following steps to avoid the direct computation of VShA. 
Firstly, we determine the solution to the following equation which we will denote p afterward 



/ 2((2 - 0)wi(|VAp) + ((/)- l)v'2{\\7A\^))Vp ■ VAVA • V^-da; 
Jn 

+ [ u(x, |VA|^)Vp- ViAdcc^ / VShA- {V A - M)dx. 
Jn Jr 



(18) 



Secondly, by letting the test function for the variable (p and tp in (fT7| and (fT8| being p and ShA, respectively, we can get 
the following equation 

- / iv2i\VA\^)-vi{\VA\^))hVA-Vi;dx^ [ VShA ■ {V A -M)dx. (19) 
Jn Jr 

By comparing Eq. (|14p with Eq. ([T5)) . we obtain 

ShFi = - / {v2{\VAf ) - vi{\VAf ))hVA ■ Vpdx. (20) 
Jn 

After deriving the formula of ShFi, we simply obtain DFi by simply projecting (t;2(|VAp) — vi{\VA\'^))VA ■ Vp onto 
the finite element space Wg^'^(ri) by solving the following equation 



' / {v2{\VA\^) ~vi{\VA\^ j)VA-Vp(pdx= [ DFiipdx, \Iip ^W^ ^in). 
Jn Jn 



(21) 
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In the numerical tests, it often comes out that the level set function can not approximate 4>exact correctly at the 
boundary of fl. To approximate the function 4>exact more accurately, we introduce a Tikhonov stabilizing term jV^pda;, 
then the constrained optimization problem (|13p becomes 

minF, subject to /!:((/))= 0, (22) 



where F = Fi -\- a J^^ |V0pda; and the coefficient a is a positive real number. 

Similar to the deduction of DFi , the gradient DF of F is calculated by solving the following equation 

/ {v2{f^A^)-vi{\VAf ))VA-Vpipdx+2a 1 Vcf) ■ Vpipdx^ [ DFipdx,Vip eWQ''^{n). (23) 
Jn Jn Jn 

In this paper, we use the Lagrangian multiplier method to convert the problem (|22[) into the following unconstrained 
optimization problem 

L{(b) = F{(b) + f h{x)K{^)dx, (24) 



where li is the Lagrangian multiplier, a /^-integrable function defined on the domain 51. 

According to the general theory of optimization, the level set function </> that we seek is the saddle point of the functional 
L((/)), that is, 

£=S+^^(2<^-3)=o, 

dL (2^) 
|^i.'(0)^O. 

By multiplying two sides of Eq. ^ by (20 — 3) and making use of the constraint (0 — 1)(0 — 2) = 0, we get the formula 
to update the multiplier 

dF 

Zi = -(20-3) — . (26) 



Substituting ([SG]) into (f25| . we have 

FiT r)F 

^ = -4(0-l)(0-2)^. (27) 



We introduce the artificial time variable t and update the level set function according to the following scheme 

{x,t) = (t)o{x) in il, 



at - d4> ™ X ±1 , ^28-) 



until the level set function </> satisfy = 0. 

To discretize (pS)) . we use the forward Euler scheme 



dL 

0fc+i = f/'A; - Atfc — 1^=0^, fc = 0, 1,2,---. (29) 

During the iteration process, for the sake of numerical stability, wc let the time step At^ satisfy the Courant-Friedrichs-Lewy 
condition 

dL 

Atfe = (TVmax| — (0fc)|, (30) 
x&n o(p 

where a € (0, 1) and h is mesh size. 

We now give the choice of 4>o and the projection used in the iteration process. We let (po be a function whose value at 
every point is neither 1 nor 2, or let (/jq be a constant function and the constant is a number between 1 and 2 (1 and 2 are 
excluded). Considering the fact that the value of the final 4> at every point x G should be either 1 or 2, wc project (j> in 
the following way: 

f 1, 0<1, 

P{i,2}W=< 2, <^>2, (31) 

I (h, otherwise, 
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after updating (j> by (f29|) at each step. 

From Eq. (P7)) . the gradient ^ is equal to when the value of is cither 1 or 2. This sometimes causes the iterative 
process unable to start or proceed, thus we exclude 1 and 2 from being the candidates of the constant for c/jq- In order to 
avoid this phenomenon in the iteration process, we count the number N of the points at which the function take 1 or 2 
after projecting the level set function (j> by Pi, 2- If the number N is equal to Nt which denotes the number of the total 
points, we stop and exit the iterative process. 

In our numerical tests, there is no such case that the iteration process stops for the reason that N = Nt- Considering 
the fact that Fi becomes smaller as the iteration proceeds, we introduce a variable osci to denote the times that the Fi 
oscillates and exit the iterative process when osci reaches a given number. 
Now we will present the piecewise constant level set algorithm (PCLSA). 

Algorithm 1 PCLSA 

Initialize 0o as a suitable function, ~ 10000 , osci = and iV = 0, Compute Nt- For fc = 0, 1, 2, ■ • ■ , 

Step 1. Use <f)k to update v : Vk{x) = vi.k{2 — (pk) + W2,fe(0fe — 1) and obtain Ak and WAk by solving Eq. (fTT|) with the 
Newton-Raphson method. 

Step 2. Compute Fi^k by p^ . let osci := osci + 1 if Fi,k > Fi,fe-i. If osci is equal to the predetermined value, exit the 

iterative process; Otherwise, go to Step 3. 
Step 3. Solve the equation to compute the gradient of Fi according to ^fe and VAk, and update §^|</)fc. 
Step 4. Set Atk by (|30|) and use the scheme (p9|) to update the level set function 0. 

Step 5. Project (j) by ((3T|) . Check the projected and obtain N. Exit the iterative process if = Nt- Or else, set 
fc := fc + 1 and 4>k = 0, go to Step 1. 

Remark 1. After exiting the iterative process, the value of the level set function at every point x ^ usually doesn't 
satisfy the constraint (0 — — 2) = 0. To make cj) satisfy the constraint, we project as below: 

after exiting the iteration process. 



5 Numerical results 

In this section, four examples are solved by Algorithm 1 (PCLSA). The workpiece we studied in the numerical tests is made 
of the hard steel and possibly contains the crack which is filled with the air. The magnetic permeability of air is close to 
1, i.e. /Lti = 1, so wi = 1. For the function V2, which is relevant to the magnetic permeability of the hard steel, is defined as 

V2{s)=di + ^——-, 



ui -r 

where the concrete values for the four variables in V2 are set the same as [5]: 

ai = 0.5, bi = 4, ci = 3, di = 0.2. 

For all examples, the domain 51 = [—0.5, 0.5] x [—0.5, 0.5] and all the numerical tests are run on the PC with Intel Core 
2 Duo 2.10 GHz processor and 2 GB RAM by the software Matlab 2010b. When solving the equation (ITT|) and ([231), we 
divide the domain fl into some rectangles with the size = hy = h = 1/dim, where dim is the number of rectangles in the 
a;— direction. For the measurement data M, we generate it in the following steps: we first find a level set function (pexact 
to represent the shape of D accurately, then solve the equation (fTTj) with the Newton-Raphson algorithm to obtain the 
solution Af,xacti finally assign VAei^act to M. In order to see the closeness of the function to ^exact, we call the built-in 
'contour' command in Matlab to plot the interface between two sub-domains after the test. For the command 'contour', we 
set the parameter for the number of the contour line as 1 and choose the red-solid-linc and the blue-dotted-line to denote 
the interface between two sub-domains represented by (jjexact and (f>, respectively. In the figures that depict the evolution 
of the level set function, the red part and the blue part represent the hard steel and the air, respectively. In the figures 
that depict the evolution of the level set function, the red part and the blue part represent the hard steel and the air. 
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respectively. In some figures, the shapes of some graphics don't look like shapes they should be. For instance, the shape 
of the circle in Example 1 is more like an octagon. In order to measure the misfit between the computed (p ^^^d 4>exact 
precisely, we count the number of the points where two functions take different values. 

In this section, we do two groups of numerical tests to show the effectiveness and robustness of PCLSA. The first group 
of tests are based on the measurement data M without noise. In order to show the flexibility of the choice of the initial 
guess of (/>, we set the function as some different initial values and observe the final shapes. The second group of tests 
are based on the measurements data M with a certain level of noise. As to the robustness of PCLSA, our algorithm can 
reconstruct the shape precisely when the noise level is up to 15%, which is superior to the algorithm in [3]. 

Example 1. In this example. We use the same configuration as in [5]. The exact shape of the domain Z? is a circle with 
the center (0.2,0.15) and the radius 0.1. And the induced current density J is defined as follows: 

Ji, V> 0.4, 

J{x,y)={ -Ji, y<-QA, (33) 
0, otherwise 

where Ji = 500. This choice of J describes the case that the workpiece is wrapped by the wires. 

For this test, we set dim = 50, a — 0.9 and a — 0.001. The initial value of the level set function is chosen as i/io = 1.5 
and the upper bound of osci is 10. The iterative process stops after 304 steps. We present the results in Fig. [T] 



In Fig. 1(a) the red solid line and the blue dotted line coincide, that is to say, PCLSA can identify the shape of the air 



gap completely. In this test, the two functions, cf) ^^nd 4>exact, take the same value at every point a; G 51. Fig. 1(b) shows 



the change of Fi with respect to the number of iterations. From Fig. 1(b) we can see that the value of Fi first becomes 
smaller, but doesn't become smaller any more after 300 steps. So it is reasonable for us to choose osci — 10. Fig. 1(c) to 



Fig. 1(f) show the evolution of the function (j) and the interface. During the iterative process, the values of (f) at the points 
of U, develop towards to our expectation except four corner points. As the increase of the number of iterations, the values 
of <j) at these four points first become small, which is away from our expectation, but becomes large after some steps. In 
Fig. [21 the values of the level set function (j) at the four corner points and sonic points at the boundary of are 1 instead 



of 2. By comparing with the picture of 4) after 300 steps (see Fig. 1(f) and Fig. [5]), we can see that it's the effects of the 
regularization to let the values of (p at four points become large. 

Example 2. In this example, the shape of the air gap in the workpiece is non-convex, which is more like a shrimp. The 
current density function J is defined as p3|) with the constant Ji equal to 500. This numerical test describes the case that 
the workpiece is all wrapped with the wires. 

In this test, we set some variables as follows: 

dim = 40, Of = 0.001, a = 0.9. 

The initial guess of (j) and the upper bound of osci are chosen as 0o = 1.5 and 15, respectively. The iterative process stops 
after 212 iterations. We present the numerical results in Fig. |3l 



In Fig. 3(a) the red solid line coincides with the blue dotted one. And the number of the points where (p (pexact 
take different values is 0. That is, the algorithm can identify the shape of D accurately. Fig. |3(b)| depicts the change of 
Fi with respect to the number of iterations. In Fig. |3(b)[ unlike Example [1] the value of Fi is below 10~^^ after some 
oscillations. Before exiting the iterative process, the times that the value of F oscillates do not exceed 15, so it's reasonable 



to choose the upper bound of osci to be 15. Fig. 3(c)| to Fig. |3(f)| describe the evolution of 4> and the interface. Similar to 



the evolution of the level set function (p in the Example [l] as the increase of the number of iterations, the values of (p at 
four corners first become smaller then start to increase to 2. In Fig. IH the values of (p at four points of O are 1. Thus we 
can conclude that it is the effects of the regularization to pull the values of (p at these points back to 2. 

Example 3. In this example, in order to test the ability of identifying the crack that is disconnected, we let the crack be 
the union of two circles and an ellipse. The current density function is chosen as J = 500. 

In this test, the variables are chosen as follows: 

dim = 50, a = 0.001, a = 0.9. 
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Figure 1: Numerical results for Example 1. (a) The final reconstruction of cp. (b) The change of the value of Fi vs. the 



Figure 2: The figure of computed without regularization. 



The initial value of the level set function and the upper bound of osci are (j)o = 1.5 and 10, separately. The iterative process 
stops after 306 steps. We present the numerical results in Fig. [5j 

In Fig. 5(a)[ the interfaces between steel and air, represented by the red solid line and the blue dotted line, coincide. 



And the number of the points which (j) and (jiexact take different values is 0. Fig. 5(b) depicts the change of Fi with respect 



to the number of iterations. In Fig. 5(b) the value of Fi first decreases and then oscillates in the last few steps. So it 



is reasonable for us to choose 10 as the upper bound of osci. Fig. 5(c) to Fig. 5(f) describe the evolution of <f> and the 
interface. As the increase of the number of iterations, the values of (j) at four corners first become smaller then start to 
increase to 2. In Fig. [6l the values of (p at four points of Vt are 1. Thus we can conclude that it is the effects of the 
regularization term to pull the values of (j) at these points to 2. 

Compared with the algorithm in [9], PCLSA takes much less steps, but reconstructs a more accurate shape. In [9], it 
takes 574 iteration steps to stop the algorithm. In our PCLSA, we only need half of that iterations, that is, totally 307 
iterations to stop the algorithm. Moreover, with more iterations in [9] the final reconstruction shape still deviates a lot 
from the exact shape. In our final result (see Fig [5]), however, there is no deviations at all. Thus PCLSA performs better 
than the algorithm in [9] for this example. 

In order to show the flexibility of the initial guess of 0, we assign different values to and compare the final results, we 
set (f)Q to be different constants, that is, 1.2, 1.4, 1.7, 1.9, respectively, and present the final results in Fig. [71 For these four 
different values of 0Oi the interfaces represented by and 4>exact coincide and the number of the points where and 4>exact 
take different value is 0. This indicates that our algorithm can identify the shape exactly. Furthermore, we choose a more 
general function as the initial value of and also observe the final interface. The general function is 0(x) = 1 + rand{x) 
where rand{x) can produce pseudo-random values between and 1. We present the numerical results in Fig. [8l where the 
interfaces represented by (j) and 4>exact coincide. 

From the above observation, PCLSA relies little on the initial value of the level set function. Thus, we need not to put 
too much attentions on the initial guess of (j) in our PCLSA for the nonlinear electromagnetism recovery problems. 

Example 4. To test the robustness of PCLSA, we use the measurement data of M with a certain level of noise. For this 
test, we use the same case in example 2 except that the measurement data of M is polluted by the certain level of noise, 
5%, 10%, 15% and 20%. 

In this test, dim = 50. We set (f>o = 1.5, the upper bound of osci to be 10, a = 0.1 and a ~ 0.9. For different noise level 
cases, the final results are shown in Fig. |9l When the noise level is under 15%, the PCLSA reconstructs the shape of D 
completely. That is, in our numerical test, the two functions and (pexact take the same value at every point a; S 51. But 
when the noise level is up to 20%, though the interface between the steel and air represented by <j) is recovered, there are 
some flaws at the interface between D and ^\D (some blue points). For the range of cr, the interval 0.6 ^ 0.9 is shown by 
our numerical tests to be not a bad choice. 

However, a affects the iteration steps in our PCLSA. The larger a is. the less steps the algorithm takes. But at high 
noise level, the reconstructed shape is better when the value of a is smaller. 
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Figure 4: The figure of witliout reguiarization. 



6 Conclusion 

In this paper, we propose a piecewise constant level set algorithm for the nonlinear electromagnetism inverse problem. In 
our numeric test, our PCLSA do not rely on the initial guess of the level set function (/), and can reconstruct the shape of 
fl exactly even when the noise level is high. Our PCSLA is quite effective and robust to solve a kind of nonlinear inverse 
problems. 
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